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We report on a new variant of the hybrid Monte Cario algorithm employing a polynomial ap- 
proximation of the inverse of the non-Hermitian Dirac-Wilson operator. Our approximation relies 
on simple and stable recurrence relations of complex Chebyshev polynomials. First performance 
figures are presented. 
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1. Introduction 

Despite the steady progress in developing new machines providing more and more computa- 
tional resources, two-flavor dynamical fermion simulations remain challenging. Moving on to large 
volume simulations we experienced frequent occurrences of large energy violations within the 
hybrid Monte Cario (HMC) update [Q] indicating possible algorithmic instabilities and raising wor- 
ries about reversibility violations. In it is pointed out that such instabilities are caused by tiny 
eigenvalues of the Dirac-Wilson operator. Therefore we propose a new HMC variant reviving the 
idea to reweight observables and approximate the inverse, non-Hermitian Dirac-Wilson operator 
by Chebyshev polynomials. Our new variant is demonstiated using Wilson's lattice action iQ] 

, V/, VA) = - £ Tr { 1 - t/p} + £ \lf{x)M^. Y{y) , (1.1) 

where the first sum runs over all plaquettes Up (x) = U^{x) ■Uv{x + fL) -U^ix + vY -Uyixy of the 
gauge field U^{x) and the second sum over all lattice sites, go is the strong coupling constant and 
M is the non-Hermitian Dirac-Wilson operator given in matrix notation. After integrating out the 
quark fields (anti-commuting Grassmann variables) ijf and y we arrive for the second summand 
at the fermion determinant det{M}. Due to the sparse structure of M it is suitable for even-odd 
preconditioning, which leads to a factorization of det{M} 

det{M} = det{Mee} • detiM"^} (1.2) 
= det{Mee} • det{M^} • det{Moo}. (1.3) 

The two possibilities are commonly named asymmetric and symmetric with the preconditioned 
operator = Mqo - Moe M¿ Meo and = 1- M¿ Moe M^ Meo, respectively. 

Setting up an algorithm to simúlate two flavor QCD, we consider the determinant of MM\ 
The algorithmic concept of our new HMC variant foUows the concept of the polynomial hybrid 
Monte Cario (PHMC) by Frezotti and Jansen [^]. They approximate the Hermitian operator by a 
root factorization and allow to compénsate by a reweighting factor for a possible deviation from 
importance sampling. Here we make use of the non-Hermitian, even-odd preconditioned operator 
M = or M^, introduce polynomials Pn ~ M^^ and incorpórate as well a reweighting factor 

áet{MM^} = det{[MP„] [MP„\^] ■ [áet{PnPl]Y\ (1.4) 

In the next section we first motívate our cholee of the non-Hermitian Dirac-Wilson opera- 
tor and show some of its properties. Our approximation of the inverse Dirac-Wilson operator in 
terms of Chebyshev polynomials is presented in Section ^ and in Section^ we introduce the basic 
steps of our Non-Hermitian Polynomial Hybrid Monte Cario (NPHMC). The dependence on the 
polynomial parameters is analyzed in Section |[ 



2. Non-Hermitian Dirac-Wilson Operator 

Using matrix notation and the hopping parameter representation we can write the Dirac-Wilson 
operator as M^y = 5xy — Kxy, where all interactions are contained in 

K-xy = K" ( ílxy — ■^'^sw'^llV'^IlvSxy ] ■ (2-1) 
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Kxy is thus built up by the hopping operator Hxy and the 0{a) improvement, the Sheikholeslami- 
Wohlert-term, and is proportional to the hopping parameter K. In we present spectral properties 
of this operator in a set-up with Schrodinger functional (SF) boundary conditions (BC). For our 
further considerations the following properties are important: 

• M has a complex spectrum lying in the positive half-plane bounded by an elhpse 

• M^^ can be approximated recursively by simple and stable polynomials 

• We expect a "good" approximation with a lower degree polynomial than if approximating 
the Hermitian operator, which is indicated by a lower condition number^ 

• These properties carry over to the even-odd preconditioned operator M 

The latter manifests itself in an exact relation between the eigenvalues X of the preconditioned 
operator K and the eigenvalues A of A' if the O {a) improvement is switched off (csw = 0) 

X{k)=X'^{K). (2.2) 

Assuming next the spectrum of K and K, respectively, to fiU perfectly an ellipse we parameterize 



the eigenvalues oí Khy X = ecosh(í? + /0) with eccentricity e and "angles" ú and 0. From ( P^ 
foUows now for the spectrum of the preconditioned operator 

é = 8 = (2.3) 

where é is the eccentricity of the ellipse bounding the spectrum of K and 5 marks a positive shift 



along the real axis. With 0{a) improvement tumed on, (2.2) as well as the drawn conclusions hold 
only approximately. 

Computing the spectral boundary numerically by using the complex Lanczos method we can 
nicely visualize the shift and the advantage due to even-odd preconditioning (see Fig. |l]). Moreover 
we see that with Sheikholeslami-Wohlert term the symmetric versión of even-odd preconditioning 
is superior because it leads to a more compact and round spectrum. In addition one can verify 
certain symmetries of the operator (see e.g. [^): without the clover term and if all dimensions of 
the lattice are even, the spectrum exhibits a mirror symmetry under sign flip and has complex pairs 
of eigenvalues (75 Hermiticity). With clover term the first is not present, but we still find complex 
pairs of eigenvalues. 

3. Chebyshev Approximation 

In order to find a polynomial approximation of the inverse, non-Hermitian Dirac-Wilson oper- 
ator M we first introduce a "small quantity" (remainder) 

Rn+,{M) = l-MPn. (3.1) 

By construction Rn+i is small on an elliptical región of the spectrum of M and suitable for an 



approximation by scaled and translated Chebyshev polynomials as introduced by Manteuffel [10] 



R.,.{M) = ^f^,- (3.2) 
Tn+]_[d/e) 



We thank Tony Kennedy for pointing that out. 
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Figure 1: Computing the spectral boundary of II ~M (upper plots) and 1 —M (lower plots) by the Lanczos 
method on a set of 50 pure-gauge configurations: 8"*^ lattice at j3 = 6.0 and K = 0.13458 in the SF. Left 
without / right with clover-term; lower right plot shows in red asymmetric and in blue symmetric even-odd 
preconditioning. 



Here d = \ + o is ideally such that the spectrum of K is origin centered and henee M = d — K. The 



spectral región approximated is bounded by an ellipse with eccentricity e. Eq. ( |3.2| ) together with 
the well-known recurrence relation for Chebyshev polynomials, 



Tn+ 1 (z) = '^zTn (z) -Tn-i (z) with Ti = z and Tq = 1 , 



(3.3) 



provides the key to obtain first a recursive description for the Rn+\ and by (3.1 ) also for our sought 



polynomials Pn [11], 



/?„_^i = ünKRn + (1 - dan)Rn-\ with 7?i = k/d and /?o = H, (3.4) 

Pn = an{'^+kPn-\) + {'\--dan)Pn-2 with Pi = ( II +^/í/) and Po = t/d, (3.5) 



The real coefficients a„ are given by 

a„ = (íí— a„_ie^/4)^' with a\ = d{d^ — /2)^^ . 



(3.6) 



and converge to lim„^oo = lid — e\j d} j — X) ¡ . Concerning the two recurrence relations 
( Pl ) and ( ^ ) we like to emphasize: although given as matrix relations, ( ^^ and (3.5) lead to 
repeated, numerically cheap matrix times vector multiplications. Furthermore, the inverting poly- 



nomial is obtained from numerically stable and simple two-step recursions.[12] 
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4. NPHMC update 

To genérate configurations with the Boltzmann weight exp{— 5} by a NPHMC update, we 
manipúlate the action ( pTT| ) accordingly 

5= — yTr{l-[/p} + 2[lndet{Mee} + lndet{Moo}] + (/>VJ(M)P„(M)(í», (4.1) 

^0 y ^ V ' 

where from now on we choose M = M since M is superior to M . The first term is the unchanged 
contribution from the gauge fields, it foUows the determinant contribution due to even-odd precon- 
ditioning and the last term gives rise to the bosonic contribution obtained by estimating the second 



factor of the right hand side in (1.4) as bosonic integral. We take care of the first factor in (1.4), 
det{[M/'„][M/',¡]^}, by calculating a stochastic estímate using random Gaussian fields r]c 



C = exp {t]¿[1I -{P¡MmPn)-']r]c} ■ 



(4.2) 



This term is only computed when measuring observables. The general procedure of an HMC 
update is unchanged. Differences occur in the way the bosonic contribution 5^, is treated. 

At the beginning of a trajectory the pseudo-fermion fields have to be generated with the 
correct distribution by inverting P„ on a random Gaussian field 77 



(II-/?, 



n+l) 



^-1 



Mrj. 



(4.3) 



Here it is important to use the second relation and by that invert a well conditioned matrix needing 
only a little number of expensive iteration steps. Next we compute the variation d{Sh} as it is 
required to intégrate the equations of motions. Varying each occurrence of K and reorganizing the 
expressions we yield 



S{Sb} = '£ ^l_iaid{K}x,-r+U.c. 



(4.4) 



with 



Xj = aj^ + ajKXj-i + {l-daj)Xj-2; Xi = ai{í+K/d)(¡); Xo = ^/'^ (4-5) 

C+J = C+j-Mn-j+l + C+;-2(l - dan-j+l); C+l = Ú^^n\ ^¡ = xl (4-6) 

Finally, the reweighting factor C is estimated requiring again the inversión of a well-conditioned 



matrix. Of course also here we can replace MPn by II — For further details see [12]. 



5. Performance Tests 



Testing the proposed new HMC variant we investígate the dependence on the three polynomial 
parameters 5, é and n on an 8^ lattice at j8 = 6.0 and K = 0.13458. The tests are performed by 
keeping two parameters fixed, while varying the third one and monitoring as observables: the mean 
valué of the correction factor (C), the relative width of its distribution qc = \l {C^) — {C)^/ (C) as 
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Figure 2: Dependence of the correction factor C and the width of its distribution C£ on 5 (light blue o) and 
é (dark blue A), respectively. 
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Figure 3: Dependence of the number of CF iterations on 5 (light blue o) and é (dark blue A), respectively. 



well as the number of conjúgate gradient (CG) iterations required to compute the correction factor, 
#(iterations CG). 

As can be seen by looking at Figs. ^ and |3| the two parameters specifying the elliptical región 
of the approximating polynomial favor e « 5 « 0.36 but do not require a sophisticated fine tuning. 
This valué is in good agreement with the ones concluded from the bounding ellipse shown in Fig. |l] 
in case of symmetric even-odd preconditioning. 

Tuming to the dependence on the polynomial degree n (Fig. ^ we find a much stronger de- 
pendence. This parameter is mainly responsible for the quality of our approximation. Henee it 
affects the cost as well as the noisiness of the correction factor. For a higher degree polynomial 
C approaches one, while qc and #(iterations CG) go to zero. Using the 4n ■ #(iterations CG) as a 
coarse estímate for the cost, a good cholee of the polynomial degree - for this set-up -isn = 50. 

6. Conclusión 

Employing the scaled and translated Chebyshev polynomials to approximate the inverse, non- 
Hermitian Dirac-Wilson operators allows to derive a variant of the HMC update algorithm based 
on simple and stable recurrence relations. The approximation does not require a fine tuning of 
the parameters specifying the elliptical approximation región, whereas the polynomial degree n is 
important for the quaUty of the approximation and the numerical cost. 
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Figure 4: Dependence of C (A), (^) ^nd #(iterations CG) (o) on the polynomial degree n. 

A conclusive comparison between the performance of different HMC-type algorithm has not 
been performed yet. However, the first data indícate that the 1-pseudo-fermion NPHMC is shghtly 
superior than a standard 1-pseudo-fermion HMC, but inferior to an HMC versión incorporating 



the Hasenbusch-trick [13] and múltiple time scale integration [14|. Constructing a NPHMC with 
two pseudo-fermions by means of the Hasenbusch-trick is possible. Unfortunately, this forces an 
involved tuning of the polynomial degrees and appears to be not too promising. 

As a byproduct of our studies of the non-Hermitian Dirac-Wilson operator, we found that the 
symmetric even-odd preconditioned operator is advantageous because of a more compact spectrum. 
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